/*!
 * \file numerics_template.cpp
 * \brief This file contains all the convective term discretization.
 * \author F. Palacios
 * \version 5.0.0 "Raven"
 *
 * SU2 Lead Developers: Dr. Francisco Palacios (Francisco.D.Palacios@boeing.com).
 *                      Dr. Thomas D. Economon (economon@stanford.edu).
 *
 * SU2 Developers: Prof. Juan J. Alonso's group at Stanford University.
 *                 Prof. Piero Colonna's group at Delft University of Technology.
 *                 Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology.
 *                 Prof. Alberto Guardone's group at Polytechnic University of Milan.
 *                 Prof. Rafael Palacios' group at Imperial College London.
 *                 Prof. Edwin van der Weide's group at the University of Twente.
 *                 Prof. Vincent Terrapon's group at the University of Liege.
 *
 * Copyright (C) 2012-2017 SU2, the open-source CFD code.
 *
 * SU2 is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * SU2 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
 */

#include "../include/numerics_structure.hpp"
#include "../include/fluid_model.hpp"
#include <limits>

double calcLmd(double aRho, double aU, double aV, double aE, std::size_t aK);
double calcg(double aU, double au, double aV, double av, double aRho, double aLmd);
double sign(double x);
double calcErfc(double x);
void calcUPsi(double aRho, double aU, double aLmd, int aUFlag, int maxIndex, double* aUPsi);
void calcUmVnXPsi(double* aUPsi, int aUm, double* aVPsi, int aVn,  
        const double* aArr, double aKsi2, double aKsi4, double* aValue);

double* calcFlux(bool isINV, CFluidModel* aFluidModel, double aGas_Constant, double aPrandtl_Lam, double* aVar, double** GradientA, double* aCG, 
        double* bVar, double** GradientB, double* bCG,
        double* sfCG, double surfaceVolume, double* aNormal, 
        int nDim, int nVar, double aK, double deltaT, double eddyVis, double* lmdij);
void solveEquation(double aPrandtl, double aRho, double aU, double aV, double aLambda, double* B, double aK, double* aX);
void transformCoordSystem(double* aCG, double* aVar, double** aGradient, 
        double* rCG, double* rVar, double** rGradient, double* normal, double* tangent);
void transformConsVar(double* aVar, double* rVar, double* normal, double* tangent);
void transformConsGrad(double** aGradient, double** rGradient, double* normal, double* tangent);
void transformCoord(double* aCoord, double* rCoord, double* normal, double* tangent);

CGasKineticScheme::CGasKineticScheme(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) {

    implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
    grid_movement = config->GetGrid_Movement();
    kappa = config->GetRoe_Kappa(); // 1 is unstable

    Gamma = config->GetGamma();
    Gamma_Minus_One = Gamma - 1.0;

    Diff_U = new su2double [nVar];
    Velocity_i = new su2double [nDim];
    Velocity_j = new su2double [nDim];
    RoeVelocity = new su2double [nDim];
    delta_vel  = new su2double [nDim];
    delta_wave = new su2double [nVar];
    ProjFlux_i = new su2double [nVar];
    ProjFlux_j = new su2double [nVar];
    Lambda = new su2double [nVar];
    Epsilon = new su2double [nVar];
    P_Tensor = new su2double* [nVar];
    invP_Tensor = new su2double* [nVar];
    for (iVar = 0; iVar < nVar; iVar++) {
        P_Tensor[iVar] = new su2double [nVar];
        invP_Tensor[iVar] = new su2double [nVar];
    }

    MeanVelocity = new su2double [nDim];
    PrimVar_i = new su2double [nDim+3];
    PrimVar_j = new su2double [nDim+3];
    Mean_PrimVar = new su2double [nDim+3];

    Proj_Mean_GradPrimVar_Edge = new su2double [nDim+1];
    Mean_GradPrimVar = new su2double* [nDim+1];
    for (iVar = 0; iVar < nDim+1; iVar++)
        Mean_GradPrimVar[iVar] = new su2double [nDim];
    Edge_Vector = new su2double [nDim];
}

CGasKineticScheme::~CGasKineticScheme(void) {
    unsigned short iVar;

    delete [] Diff_U;
    delete [] Velocity_i;
    delete [] Velocity_j;
    delete [] RoeVelocity;
    delete [] delta_vel;
    delete [] delta_wave;
    delete [] ProjFlux_i;
    delete [] ProjFlux_j;
    delete [] Lambda;
    delete [] Epsilon;
    for (iVar = 0; iVar < nVar; iVar++) {
        delete [] P_Tensor[iVar];
        delete [] invP_Tensor[iVar];
    }
    delete [] P_Tensor;
    delete [] invP_Tensor;

    delete [] MeanVelocity;
    delete [] PrimVar_i;
    delete [] PrimVar_j;
    delete [] Mean_PrimVar;
    delete [] Proj_Mean_GradPrimVar_Edge;
    delete [] Edge_Vector;
    for (iVar = 0; iVar < nDim+1; iVar++)
        delete [] Mean_GradPrimVar[iVar];
    delete [] Mean_GradPrimVar;
}

void CGasKineticScheme::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) 
{
    su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0};
    su2double ProjGridVel = 0.0;


    /*--- Pressure, density, enthalpy, energy, and velocity at points i and j ---*/

    Pressure_i = V_i[nDim+1];                       Pressure_j = V_j[nDim+1];
    Density_i = V_i[nDim+2];                        Density_j = V_j[nDim+2];
    Enthalpy_i = V_i[nDim+3];                       Enthalpy_j = V_j[nDim+3];
    SoundSpeed_i = V_i[nDim+4];                     SoundSpeed_j = V_j[nDim+4];
    Energy_i = Enthalpy_i - Pressure_i/Density_i;   Energy_j = Enthalpy_j - Pressure_j/Density_j;

    sq_vel_i = 0.0; sq_vel_j = 0.0;
    for (iDim = 0; iDim < nDim; iDim++) {
        Velocity_i[iDim] = V_i[iDim+1];
        Velocity_j[iDim] = V_j[iDim+1];
        sq_vel_i += 0.5*Velocity_i[iDim]*Velocity_i[iDim];
        sq_vel_j += 0.5*Velocity_j[iDim]*Velocity_j[iDim];
    }

    /*--- Recompute conservative variables ---*/

    U_i[0] = Density_i; U_j[0] = Density_j;
    for (iDim = 0; iDim < nDim; iDim++) {
        U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim];
    }
    U_i[nDim+1] = Density_i*Energy_i; U_j[nDim+1] = Density_j*Energy_j;

    Area = 0;
    for (iDim = 0; iDim < nDim; iDim++)
        /*!< \brief Normal: Normal vector, it norm is the area of the face. */
        Area += Normal[iDim]*Normal[iDim];
    Area = sqrt(Area);                    /*! Area of the face*/

    for (iDim = 0; iDim < nDim; iDim++)
        UnitNormal[iDim] = Normal[iDim]/Area;   /* ! Unit Normal*/

    /*--- Compute mean values of the variables ---*/

    R = sqrt(fabs(Density_j/Density_i));
    RoeDensity = R*Density_i;

    sq_vel = 0.0;
    for (iDim = 0; iDim < nDim; iDim++) {
        RoeVelocity[iDim] = (R*Velocity_j[iDim]+Velocity_i[iDim])/(R+1);
        sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim];
    }
    RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1);

    RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel);

    /*--- Negative RoeSoundSpeed2, the jump 
      variables is too large, exit the subrotuine 
      without computing the fluxes ---*/

    if (RoeSoundSpeed2 <= 0.0) {
        for (iVar = 0; iVar < nVar; iVar++) {
            val_residual[iVar] = 0.0;
            for (jVar = 0; jVar < nVar; jVar++) {
                val_Jacobian_i[iVar][iVar] = 0.0;
                val_Jacobian_j[iVar][iVar] = 0.0;
            }
        }
        return;
    }

    RoeSoundSpeed = sqrt(RoeSoundSpeed2);

    /*--- Compute P and Lambda (do it with the Normal) ---*/

    GetPMatrix(&RoeDensity, RoeVelocity, &RoeSoundSpeed, UnitNormal, P_Tensor);

    /*--- Compute ProjFlux_i ---*/

    GetInviscidProjFlux(&Density_i, Velocity_i, &Pressure_i, &Enthalpy_i, Normal, ProjFlux_i);

    /*--- Compute ProjFlux_j ---*/

    GetInviscidProjFlux(&Density_j, Velocity_j, &Pressure_j, &Enthalpy_j, Normal, ProjFlux_j);

    ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0;
    for (iDim = 0; iDim < nDim; iDim++) {
        ProjVelocity   += RoeVelocity[iDim]*UnitNormal[iDim];
        ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim];
        ProjVelocity_j += Velocity_j[iDim]*UnitNormal[iDim];
    }

    /*--- Projected velocity adjustment due to mesh motion ---*/

    if (grid_movement) {
        ProjGridVel = 0.0;
        for (iDim = 0; iDim < nDim; iDim++) {
            ProjGridVel   += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*UnitNormal[iDim];
        }
        ProjVelocity   -= ProjGridVel;
        ProjVelocity_i -= ProjGridVel;
        ProjVelocity_j -= ProjGridVel;
    }

    /*--- Flow eigenvalues and entropy correctors ---*/

    for (iDim = 0; iDim < nDim; iDim++)
        Lambda[iDim] = ProjVelocity;

    Lambda[nVar-2] = ProjVelocity + RoeSoundSpeed;
    Lambda[nVar-1] = ProjVelocity - RoeSoundSpeed;

    /*--- Compute absolute value with Mavriplis' entropy correction ---*/

    MaxLambda = fabs(ProjVelocity) + RoeSoundSpeed;
    Delta = config->GetEntropyFix_Coeff();



    double eddyVis = 0.0;

    unsigned short turb_model = config->GetKind_Turb_Model();
    if(turb_model) 
    {
        //cout << "flux" << endl;
        //getchar();
        double eddyVisa = V_i[nDim+6];
        double eddyVisb = V_j[nDim+6];
        /*
           double** PrimGradientA = node[aPoint]->GetGradient_Primitive(); 
           double** PrimGradientB = node[bPoint]->GetGradient_Primitive(); 
           eddyVisa = eddyVisa + PrimGradientA[nDim + 6][0] * (edgeCG[0] - aCG[0])
           + PrimGradientA[nDim + 6][1] * (edgeCG[1] - aCG[1]);
           eddyVisb = eddyVisb + PrimGradientB[nDim + 6][0] * (edgeCG[0] - bCG[0])
           + PrimGradientB[nDim + 6][1] * (edgeCG[1] - bCG[1]);
           */
        eddyVis = 0.5 * (eddyVisa + eddyVisb);

        //cout << eddyVisa << "  " <<eddyVisb << endl;
        //cout << eddyVis << endl;
        //cout << "flux" << endl;
        //getchar();

    }
    double tempTimeStep = DeltaTimeStep;

    double templmd = 0.0;
    bool isINV = config->GetGKS_NS_INV();
    double* tempflux = calcFlux(isINV, FluidModel, Gas_Constant, Prandtl_Lam, U_i, ConsVar_Grad_i, Coord_i, 
            U_j, ConsVar_Grad_j, Coord_j,
            Coord_ij, Area, UnitNormal, nDim, nVar, K, tempTimeStep, eddyVis, &templmd);

    for (iVar = 0; iVar < nVar; iVar++) {
        //Lambda[iVar] = max(fabs(Lambda[iVar]), Delta*MaxLambda);
        Lambda[iVar] = templmd;
    }

    for(iVar = 0; iVar < nVar; iVar++) val_residual[iVar] = tempflux[iVar] / tempTimeStep;

    if (implicit)
    {

        /*--- Compute inverse P ---*/

        GetPMatrix_inv(&RoeDensity, RoeVelocity, &RoeSoundSpeed, UnitNormal, invP_Tensor);

        /*--- Jacobians of the inviscid flux, scaled by
          kappa because val_resconv ~ kappa*(fc_i+fc_j)*Normal ---*/

        GetInviscidProjJac(Velocity_i, &Energy_i, Normal, kappa, val_Jacobian_i);
        GetInviscidProjJac(Velocity_j, &Energy_j, Normal, kappa, val_Jacobian_j);

        /*--- Diference variables iPoint and jPoint ---*/

        for (iVar = 0; iVar < nVar; iVar++)
            Diff_U[iVar] = U_j[iVar]-U_i[iVar];

        /*--- Roe's Flux approximation ---*/

        for (iVar = 0; iVar < nVar; iVar++) {

            //val_residual[iVar] = kappa*(ProjFlux_i[iVar]+ProjFlux_j[iVar]);
            for (jVar = 0; jVar < nVar; jVar++) {
                Proj_ModJac_Tensor_ij = 0.0;

                /*--- Compute |Proj_ModJac_Tensor| = P x |Lambda| x inverse P ---*/

                for (kVar = 0; kVar < nVar; kVar++)
                    Proj_ModJac_Tensor_ij += P_Tensor[iVar][kVar]*Lambda[kVar]*invP_Tensor[kVar][jVar];

                //val_residual[iVar] -= (1.0-kappa)*Proj_ModJac_Tensor_ij*Diff_U[jVar]*Area;
                val_Jacobian_i[iVar][jVar] += (1.0-kappa)*Proj_ModJac_Tensor_ij*Area;
                val_Jacobian_j[iVar][jVar] -= (1.0-kappa)*Proj_ModJac_Tensor_ij*Area;

            }

        }

        /*--- Jacobian contributions due to grid motion ---*/

        if (grid_movement) {
            ProjVelocity = 0.0;
            for (iDim = 0; iDim < nDim; iDim++)
                ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
            for (iVar = 0; iVar < nVar; iVar++) {
                //val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);

                /*--- Implicit terms ---*/

                val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity;
                val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity;
            }
        }

        if(isCalcVisJacobian)
        {
            for (iVar = 0; iVar < nDim+3; iVar++) {
                PrimVar_i[iVar] = V_i[iVar];
                PrimVar_j[iVar] = V_j[iVar];
                Mean_PrimVar[iVar] = 0.5*(PrimVar_i[iVar]+PrimVar_j[iVar]);
            }

            /*--- Compute vector going from iPoint to jPoint ---*/

            dist_ij_2 = 0.0;
            for (iDim = 0; iDim < nDim; iDim++) {
                Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim];
                dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim];
            }

            /*--- Laminar and Eddy viscosity ---*/

            Laminar_Viscosity_i = V_i[nDim+5]; Laminar_Viscosity_j = V_j[nDim+5];
            Eddy_Viscosity_i = V_i[nDim+6]; Eddy_Viscosity_j = V_j[nDim+6];

            /*--- Mean Viscosities and turbulent kinetic energy ---*/

            Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j);
            Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j);
            Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j);

            /*--- Projection of the mean gradient in the direction of the edge ---*/

            for (iVar = 0; iVar < nDim+1; iVar++) {
                Proj_Mean_GradPrimVar_Edge[iVar] = 0.0;

                for (iDim = 0; iDim < nDim; iDim++) {
                    Mean_GradPrimVar[iVar][iDim] = 0.5*(PrimVar_Grad_i[iVar][iDim] + PrimVar_Grad_j[iVar][iDim]);
                    Proj_Mean_GradPrimVar_Edge[iVar] += Mean_GradPrimVar[iVar][iDim]*Edge_Vector[iDim];
                }

                if (dist_ij_2 != 0.0) {
                    for (iDim = 0; iDim < nDim; iDim++) {
                        Mean_GradPrimVar[iVar][iDim] -= (Proj_Mean_GradPrimVar_Edge[iVar] -
                                (PrimVar_j[iVar]-PrimVar_i[iVar]))*Edge_Vector[iDim] / dist_ij_2;
                    }
                }

            }

            /*--- Get projected flux tensor ---*/

            GetViscousProjFlux(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, Normal, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);

            /*--- Compute the implicit part ---*/


            if (dist_ij_2 != 0.0) {

                double** temp_jacobian_i;
                double** temp_jacobian_j;

                temp_jacobian_i = new su2double* [nVar];
                temp_jacobian_j = new su2double* [nVar];
                for (iVar = 0; iVar < nVar; iVar++) {
                    temp_jacobian_i[iVar] = new su2double [nVar];
                    temp_jacobian_j[iVar] = new su2double [nVar];
                }

                GetViscousProjJacs(Mean_PrimVar, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,
                        sqrt(dist_ij_2), UnitNormal, Area, Proj_Flux_Tensor, temp_jacobian_i, temp_jacobian_j);
                for(iVar = 0; iVar < nVar; iVar++)
                {
                    for(int ii = 0; ii < nVar; ii++)
                    {
                        val_Jacobian_i[iVar][ii] -= temp_jacobian_i[iVar][ii];
                        val_Jacobian_j[iVar][ii] -= temp_jacobian_j[iVar][ii];
                    }
                }

                for(iVar=0; iVar<nVar; iVar++)
                {
                    delete[] temp_jacobian_i[iVar];
                    delete[] temp_jacobian_j[iVar];
                }
                delete[] temp_jacobian_i;
                delete[] temp_jacobian_j;

            }

        }

    }

    delete[] tempflux;

}


double calcLmd(double aRho, double aU, double aV, double aE, std::size_t aK)
{
    double tempLmd = (aK + 2.0) * aRho / (4.0 * (aE - 0.5 * aRho * (aU * aU + aV * aV)));
    return tempLmd;
};

double calcg(double aU, double au, double aV, double av, double aRho, double aLmd)
{
    return aRho * (aLmd / PI_NUMBER) 
        * exp(-aLmd * ((au-aU)*(au-aU) + (av-aV)*(av-aV)));
};

double sign(double x)
{
    double temp = 0.0;
    if (x > 0.0)
        temp = 1.0;
    else if (x < 0.0)
        temp = -1.0;
    else 
        temp = 0.0;
    return temp;
};

/*
 * calcErfc
 */
double calcErfc(double x)
{
    double sb, sc, erf;
    sc=fabs(x);
    sb=1.0/(1.0+0.3275911*sc);
    erf=1.0-sb*(0.254829592+sb*(-0.284496736+sb*(1.421413741
                    +sb*(-1.453152027+sb*1.061405429))))*exp(-sc*sc);

    if (x < 0.0)
    {
        erf *= -1;
    }

    return 1.0-erf;
} // calcErfc
/*
 * calcUPsi
 * aUFlag = 0:U;
 * aUFlag = 1:U>0;
 * aUFlag = 2:U<0;
 */
void calcUPsi(double aRho, double aU, double aLmd, int aUFlag, int maxIndex, double* aUPsi)
{
    if (aUFlag == 0)
    {
        aUPsi[0] = 1.0;
        aUPsi[1] = aU;
    } else if (aUFlag == 1)
    {
        aUPsi[0] = 0.5 * calcErfc(-sqrt(aLmd) * aU);
        aUPsi[1] = aU * aUPsi[0] 
            + 0.5 * exp(-aLmd*aU*aU) / (sqrt(aLmd*PI_NUMBER));
    } else if (aUFlag == 2)
    {
        aUPsi[0] = 0.5 * calcErfc(sqrt(aLmd) * aU);
        aUPsi[1] = aU * aUPsi[0] 
            - 0.5 * exp(-aLmd*aU*aU) / (sqrt(aLmd*PI_NUMBER));
    }

    for(int i=2; i<maxIndex; i++)
    {
        aUPsi[i] = aU * aUPsi[i-1] 
            + (i-1.0) * aUPsi[i-2] / (2.0 * aLmd);
    }

} //: calcUPsi

void calcUmVnXPsi(double* aUPsi, int aUm, double* aVPsi, int aVn,  
        const double* aArr, double aKsi2, double aKsi4, double* aValue)
{
    for(int ivar = 0; ivar < 4; ivar++) aValue[ivar] = 0.0;
    // calculating aValue[ivar]:  0 <= ivar < 3
    for(int ivar = 0; ivar < 3; ivar++)
    {
        if (ivar == 1) 
        {
            aUm++;
        } else if (ivar == 2) 
        {
            aUm--; aVn++;
        }
        aValue[ivar] = aArr[0] * aUPsi[aUm] * aVPsi[aVn]
            + aArr[1] * aUPsi[aUm+1] * aVPsi[aVn]
            + aArr[2] * aUPsi[aUm] * aVPsi[aVn+1]
            + 0.5 * aArr[3] * ( aUPsi[2+aUm] * aVPsi[aVn]
                    + aUPsi[aUm] * aVPsi[aVn+2] 
                    + aUPsi[aUm] * aVPsi[aVn] * aKsi2
                    );
    }

    aVn--;
    // calculating aValue[3]
    for(int i = 0; i < 2; i++)
    {
        if (i == 0)
        {
            // aValue[3] : 0.5 * u^2
            aUm += 2;
        } else if(i == 1)
        {
            // aValue[3] : 0.5 * v^2
            aUm -= 2; aVn += 2;
        }
        aValue[3] += 0.5 * (aArr[0] * aUPsi[aUm] * aVPsi[aVn]
                + aArr[1] * aUPsi[aUm+1] * aVPsi[aVn]
                + aArr[2] * aUPsi[aUm] * aVPsi[aVn+1]
                + 0.5 * aArr[3] * ( aUPsi[2+aUm] * aVPsi[aVn]
                    + aUPsi[aUm] * aVPsi[aVn+2] 
                    + aUPsi[aUm] * aVPsi[aVn] * aKsi2
                    ));
    }

    // aValue[3] : ksi2
    aVn -= 2;
    aValue[3] += 0.5 * aKsi2 * (aArr[0] * aUPsi[aUm] * aVPsi[aVn]
            + aArr[1] * aUPsi[aUm+1] * aVPsi[aVn]
            + aArr[2] * aUPsi[aUm] * aVPsi[aVn+1]
            + 0.5 * aArr[3] * ( aUPsi[2+aUm] * aVPsi[aVn]
                + aUPsi[aUm] * aVPsi[aVn+2] 
                ));

    // aValue[3] : ksi4
    aValue[3] += 0.25 * aArr[3] * aUPsi[aUm] * aVPsi[aVn] * aKsi4;
}


double* calcFlux(bool isINV, CFluidModel* aFluidModel, double aGas_Constant, double aPrandtl_Lam, double* aVar, double** GradientA, double* aCG, 
        double* bVar, double** GradientB, double* bCG,
        double* sfCG, double surfaceVolume, double* aNormal, 
        int nDim, int nVar, double aK, double deltaT, double eddyVis, double* lmdij)
{

    double qs = 0.0;
    int nParent = 2;
    double surfaceVar[nParent][nVar]; double surfaceLmd[nParent];
    double var0[nVar]; double nodeVar[2][nVar]; double lmd0 = 0.0;
    double surfaceCG[2];
    surfaceCG[0] = sfCG[0];
    surfaceCG[1] = sfCG[1];

    double*** GradientArr;
    GradientArr= new double**[nParent];
    for(int iParent = 0; iParent < nParent; iParent++)
    {
        GradientArr[iParent] = new double*[nVar];
        for(int ivar = 0; ivar < nVar; ivar++)
        {
            GradientArr[iParent][ivar] = new double[nDim];

        }
    }

    double a[nParent][2][nVar]; double A[nParent][nVar];
    double aBar[3][nVar]; double ABar[nVar];

    surfaceLmd[0] = 0.0; surfaceLmd[1] = 0.0;

    for(int ivar = 0; ivar < nVar; ivar++)
    {
        GradientArr[0][ivar][0] = GradientA[ivar][0];
        GradientArr[0][ivar][1] = GradientA[ivar][1];

        GradientArr[1][ivar][0] = GradientB[ivar][0];
        GradientArr[1][ivar][1] = GradientB[ivar][1];


        nodeVar[0][ivar] = aVar[ivar]; nodeVar[1][ivar] = bVar[ivar];

        surfaceVar[0][ivar] = 0.0; surfaceVar[1][ivar] = 0.0;

        var0[ivar] = 0.0;


        A[0][ivar] = 0.0; A[1][ivar] = 0.0;
        aBar[0][ivar] = 0.0; aBar[1][ivar] = 0.0; ABar[ivar] = 0.0;
    }

    double nodeCG[2][2];
    nodeCG[0][0] = aCG[0]; nodeCG[0][1] = aCG[1]; nodeCG[1][0] = bCG[0]; nodeCG[1][1] = bCG[1];

    /************************************************************************/
    // calculating the varR / varL
    //
    for(int iParent = 0; iParent < nParent; iParent++)
    {
        double* vectori = new double[nDim];
        for(int iDim = 0; iDim < nDim; iDim++)
        {
            vectori[iDim] = surfaceCG[iDim] - nodeCG[iParent][iDim];
        }

        // calculating the var[iParent]
        for(int ivar = 0; ivar < nVar; ivar++)
        {
            double project = 0.0;
            for(int iDim = 0; iDim < nDim; iDim++)
            {
                project += vectori[iDim] * GradientArr[iParent][ivar][iDim];
            }

            surfaceVar[iParent][ivar] = nodeVar[iParent][ivar] + project;
        }

        delete[] vectori;


        // calculating the lmd[iParent]
        surfaceLmd[iParent] = calcLmd(surfaceVar[iParent][0], 
                surfaceVar[iParent][1] / surfaceVar[iParent][0], 
                surfaceVar[iParent][2] / surfaceVar[iParent][0], 
                surfaceVar[iParent][3], aK);
        //cout << surfaceLmd[iParent] << " ";
    }


    if(surfaceLmd[0] <= 0.0 || surfaceLmd[1] <= 0.0)
    {
        for(int ivar = 0; ivar < nVar; ivar++)
        {
            GradientArr[0][ivar][0] = 0.0;
            GradientArr[0][ivar][1] = 0.0;

            GradientArr[1][ivar][0] = 0.0;
            GradientArr[1][ivar][1] = 0.0;
        }

        /************************************************************************/
        // calculating the varR / varL
        //
        for(int iParent = 0; iParent < nParent; iParent++)
        {
            double* vectori = new double[nDim];
            for(int iDim = 0; iDim < nDim; iDim++)
            {
                vectori[iDim] = surfaceCG[iDim] - nodeCG[iParent][iDim];
            }

            // calculating the var[iParent]
            for(int ivar = 0; ivar < nVar; ivar++)
            {
                double project = 0.0;
                for(int iDim = 0; iDim < nDim; iDim++)
                {
                    project += vectori[iDim] * GradientArr[iParent][ivar][iDim];
                }

                surfaceVar[iParent][ivar] = nodeVar[iParent][ivar] + project;
            }

            delete[] vectori;

            // calculating the lmd[iParent]
            surfaceLmd[iParent] = calcLmd(surfaceVar[iParent][0], 
                    surfaceVar[iParent][1] / surfaceVar[iParent][0], 
                    surfaceVar[iParent][2] / surfaceVar[iParent][0], 
                    surfaceVar[iParent][3], aK);
            //cout << surfaceLmd[iParent] << " ";
        }

    }

    //cout << endl;
    /*
     *transform the coordinate system
     */
    double temp1 = 0.0; double temp2 = 0.0;
    for(size_t i = 0; i < 2; i++)
    {
        // macro-variables
        temp1 = nodeVar[i][1]; temp2 = nodeVar[i][2];
        nodeVar[i][1] = temp1 * aNormal[0] + temp2 * aNormal[1];
        nodeVar[i][2] = -temp1 * aNormal[1] + temp2 * aNormal[0];

        // macro-variables
        temp1 = surfaceVar[i][1]; temp2 = surfaceVar[i][2];
        surfaceVar[i][1] = temp1 * aNormal[0] + temp2 * aNormal[1];
        surfaceVar[i][2] = -temp1 * aNormal[1] + temp2 * aNormal[0];

        // nodeCG
        temp1 = nodeCG[i][0]; temp2 = nodeCG[i][1];
        nodeCG[i][0] = temp1 * aNormal[0] + temp2 * aNormal[1];
        nodeCG[i][1] = -temp1 * aNormal[1] + temp2 * aNormal[0];
    }

    // surfaceCG
    temp1 = surfaceCG[0]; temp2 = surfaceCG[1];
    surfaceCG[0] = temp1 * aNormal[0] + temp2 * aNormal[1];
    surfaceCG[1] = -temp1 * aNormal[1] + temp2 * aNormal[0];

    // Gradient
    for(int iParent = 0; iParent < nParent; iParent++)
    {
        temp1 = GradientArr[iParent][0][0];
        temp2 = GradientArr[iParent][0][1];
        GradientArr[iParent][0][0] = temp1 * aNormal[0] + temp2 * aNormal[1];
        GradientArr[iParent][0][1] = -temp1 * aNormal[1] + temp2 * aNormal[0];

        temp1 = GradientArr[iParent][3][0];
        temp2 = GradientArr[iParent][3][1];
        GradientArr[iParent][3][0] = temp1 * aNormal[0] + temp2 * aNormal[1];
        GradientArr[iParent][3][1] = -temp1 * aNormal[1] + temp2 * aNormal[0];


        temp1 = GradientArr[iParent][1][0];
        temp2 = GradientArr[iParent][2][0];
        double temp3 = GradientArr[iParent][1][1];
        double temp4 = GradientArr[iParent][2][1];

        GradientArr[iParent][1][0] = aNormal[0] * aNormal[0] * temp1 
            + aNormal[0] * aNormal[1] * temp2
            + aNormal[0] * aNormal[1] * temp3
            + aNormal[1] * aNormal[1] * temp4;

        GradientArr[iParent][2][0] = -aNormal[0] * aNormal[1] * temp1 
            + aNormal[0] * aNormal[0] * temp2
            - aNormal[1] * aNormal[1] * temp3
            + aNormal[0] * aNormal[1] * temp4;;

        GradientArr[iParent][1][1] = -aNormal[0] * aNormal[1] * temp1 
            - aNormal[1] * aNormal[1] * temp2
            + aNormal[0] * aNormal[0] * temp3
            + aNormal[1] * aNormal[0] * temp4;

        GradientArr[iParent][2][1] = aNormal[1] * aNormal[1] * temp1 
            - aNormal[1] * aNormal[0] * temp2
            - aNormal[0] * aNormal[1] * temp3
            + aNormal[0] * aNormal[0] * temp4;
    }

    for(int iParent = 0; iParent < nParent; iParent++)
    {
        for(int iVar = 0; iVar < nVar; iVar++)
        {
            nodeVar[iParent][iVar] = nodeVar[iParent][iVar] - GradientArr[iParent][iVar][1] * (nodeCG[iParent][1] - surfaceCG[1]);
        }
        nodeCG[iParent][1] = surfaceCG[1];
    }


    /************************************************************************/
    // calc var0
    size_t maxIndex = 7;
    double UPsi[2][maxIndex];
    double VPsi[2][maxIndex];
    double Ksi2[2];
    double Ksi4[2];

    double rRho0;

    //  gl
    calcUPsi(surfaceVar[0][0], surfaceVar[0][1] / surfaceVar[0][0], surfaceLmd[0], 1, maxIndex, UPsi[0]);
    calcUPsi(surfaceVar[0][0], surfaceVar[0][2] / surfaceVar[0][0], surfaceLmd[0], 0, maxIndex, VPsi[0]);
    Ksi2[0] = 0.5 * aK / surfaceLmd[0];
    Ksi4[0] = 0.25 * (aK * aK + 2.0 * aK) / (surfaceLmd[0] * surfaceLmd[0]);

    //  gr
    calcUPsi(surfaceVar[1][0], surfaceVar[1][1] / surfaceVar[1][0], surfaceLmd[1], 2, maxIndex, UPsi[1]);
    calcUPsi(surfaceVar[1][0], surfaceVar[1][2] / surfaceVar[1][0], surfaceLmd[1], 0, maxIndex, VPsi[1]);
    Ksi2[1] = 0.5 * aK / surfaceLmd[1];
    Ksi4[1] = 0.25 * (aK * aK + 2.0 * aK) / (surfaceLmd[1] * surfaceLmd[1]);

    var0[0] = UPsi[0][0] * surfaceVar[0][0] + UPsi[1][0] * surfaceVar[1][0];
    rRho0 = 1.0 / var0[0];
    var0[1] = UPsi[0][1] * surfaceVar[0][0] + UPsi[1][1] * surfaceVar[1][0];
    var0[2] = UPsi[0][0] * VPsi[0][1] * surfaceVar[0][0] + UPsi[1][0] * VPsi[1][1] * surfaceVar[1][0];
    var0[3] = 0.5 * ((UPsi[0][2] + UPsi[0][0] * VPsi[0][2] + UPsi[0][0] * Ksi2[0]) * surfaceVar[0][0]
            + (UPsi[1][2] + UPsi[1][0] * VPsi[1][2] + UPsi[1][0] * Ksi2[1]) * surfaceVar[1][0]);

    lmd0 = calcLmd(var0[0], var0[1] * rRho0, var0[2] * rRho0, var0[3], aK);


    double rLmd0 = 1.0 / lmd0;

    double U0Psi[3][maxIndex];
    double V0Psi[maxIndex];

    // g0
    // u
    calcUPsi(var0[0], var0[1] * rRho0, lmd0, 0, maxIndex, U0Psi[0]);
    // u > 0
    calcUPsi(var0[0], var0[1] * rRho0, lmd0, 1, maxIndex, U0Psi[1]);
    // u < 0
    calcUPsi(var0[0], var0[1] * rRho0, lmd0, 2, maxIndex, U0Psi[2]);
    // v
    calcUPsi(var0[0], var0[2] * rRho0, lmd0, 0, maxIndex, V0Psi);


    double Ksi02 = 0.5 * aK * rLmd0;
    double Ksi04 = 0.25 * (aK * aK + 2.0 * aK) * (rLmd0 * rLmd0);


    //****************************************************************************************************
    // calculating the a / A
    double UPsi0[2][maxIndex];
    calcUPsi(surfaceVar[0][0], surfaceVar[0][1] / surfaceVar[0][0], surfaceLmd[0], 0, maxIndex, UPsi0[0]);
    calcUPsi(surfaceVar[1][0], surfaceVar[1][1] / surfaceVar[1][0], surfaceLmd[1], 0, maxIndex, UPsi0[1]);

    double B[4];
    double tempB[4];
    for(int iParent = 0; iParent < nParent; iParent++)
    {
        double aRho = surfaceVar[iParent][0];
        double raRho =  1.0 / aRho;
        double aU = surfaceVar[iParent][1] * raRho;
        double aV = surfaceVar[iParent][2] * raRho;

        double aLmd = surfaceLmd[iParent];

        double B[4];
        // calculating the a[iParent][0]
        for(int ivar = 0; ivar < nVar; ivar++) B[ivar] = GradientArr[iParent][ivar][0] * raRho;
        solveEquation(aPrandtl_Lam, aRho, aU, aV, aLmd, B, aK, a[iParent][0]);

        // calculating the a[iParent][1]
        for(int ivar = 0; ivar < nVar; ivar++) B[ivar] = GradientArr[iParent][ivar][1] * raRho;
        solveEquation(aPrandtl_Lam, aRho, aU, aV, aLmd, B, aK, a[iParent][1]);

        //
        // calculating the A[iParent]
        calcUmVnXPsi(UPsi0[iParent], 1,  VPsi[iParent], 0,  a[iParent][0], Ksi2[iParent], Ksi4[iParent], B);
        calcUmVnXPsi(UPsi0[iParent], 0,  VPsi[iParent], 1,  a[iParent][1], Ksi2[iParent], Ksi4[iParent], tempB);
        for(int ivar = 0; ivar < nVar; ivar++) B[ivar] = -1.0 * (B[ivar] + tempB[ivar]);
        solveEquation(aPrandtl_Lam, aRho, aU, aV, aLmd, B, aK, A[iParent]);

    }


    /*
       cout << "a[0]" << endl;
       for(size_t ivar = 0 ;ivar < 4; ivar++)
       {
       cout << setiosflags(ios::scientific) << setprecision(20) << a[0][0][ivar] << "   "
       << setprecision(20) << a[0][1][ivar]   << endl;
       }

       cout << "a[1]" << endl;
       for(size_t ivar = 0 ;ivar < 4; ivar++)
       {
       cout << setiosflags(ios::scientific) << setprecision(20) << a[1][0][ivar] 
       << setprecision(20) << "   "<< a[1][1][ivar]   << endl;;
       }

       cout << "A" << endl;
       for(size_t ivar = 0 ;ivar < 4; ivar++)
       {
       cout << setiosflags(ios::scientific) << A[0][ivar] << "   " << A[1][ivar]<< endl;;
       }
       */


    //***********************************************************************************************

    aFluidModel->SetTDState_rhoT (var0[0], 0.5 * rLmd0 / aGas_Constant); 
    double mu0 = aFluidModel->GetLaminarViscosity();

    double tau = 0.0;
    if(isINV)
    {
        tau = 0.01 * deltaT
            + fabs((surfaceVar[0][0]/ surfaceLmd[0]- surfaceVar[1][0]/ surfaceLmd[1]) 
                    / (surfaceVar[0][0]/ surfaceLmd[0] + surfaceVar[1][0]/ surfaceLmd[1])) * deltaT;
    } else
    {
        tau = 2.0 * lmd0 * rRho0 * (mu0 + eddyVis)
            + fabs((surfaceVar[0][0]/ surfaceLmd[0]- surfaceVar[1][0]/ surfaceLmd[1]) 
                    / (surfaceVar[0][0]/ surfaceLmd[0] + surfaceVar[1][0]/ surfaceLmd[1])) * deltaT;
    }

    double dist = sqrt((nodeCG[0][0] - nodeCG[1][0]) * (nodeCG[0][0] - nodeCG[1][0]) 
            + (nodeCG[0][1] - nodeCG[1][1]) * (nodeCG[0][1] - nodeCG[1][1]));
    *lmdij = fabs(var0[1]) / var0[0] + aFluidModel->GetSoundSpeed() 
        + 2.0 * (mu0 + eddyVis) / (var0[0] * dist);

    /*--------------------------------------------*/
    if(var0[0] <= 0.0 || var0[3] <= 0.0 || lmd0 <= 0.0)
    {
        cout << "\n !!! Error: SU2 has diverged. Now exiting... !!! \n" << endl;
        cout << "------------------------CGKSSolver::calcFlux2D---------------------------\n" << endl;
        cout << "var0 has negative value! \n" << endl;

        cout << "CG:"<< endl;
        for(int ii = 0; ii < 2; ii++)
        {
            cout << aCG[ii] << "  " << sfCG[ii] << "   "<< bCG[ii] << endl;
        }

        cout << "tCG:"<< endl;
        for(int ii = 0; ii < 2; ii++)
        {
            cout << nodeCG[0][ii] << "  " << surfaceCG[ii] << "  "<< nodeCG[1][ii] << endl;
        }

        cout << "normal" << endl;
        cout << aNormal[0] << "  " << aNormal[1] << endl;

        cout << "aVar, bVar:"<< endl;
        for(int iivar = 0; iivar < nVar; iivar++)
        {
            cout << aVar[iivar] << "  "<< bVar[iivar] << endl;
        }

        cout << "nodeVar0:"<< endl;
        for(int iivar = 0; iivar < nVar; iivar++)
        {
            cout << nodeVar[0][iivar] << "  "<< nodeVar[1][iivar] << endl;
        }
        cout << "Gradient" << endl;
        for(int iivar = 0; iivar < nVar; iivar++)
        {
            cout << GradientA[iivar][0] << "  " << GradientA[iivar][1] << "  "
                << GradientB[iivar][0] << "  " << GradientB[iivar][1] << endl;
        }

        cout << "var0:"<< endl;
        for(int iivar = 0; iivar < nVar; iivar++)
        {
            cout << surfaceVar[0][iivar] << "  "<< var0[iivar] << "  " << surfaceVar[1][iivar] << endl;
        }

        cout << "tau = " << tau  << "  deltaT = " << deltaT << endl;
        cout << "mu0 = " << mu0<< " eddyVis = " << eddyVis<< endl;

        exit(EXIT_FAILURE);
    }

    for(int ivar = 0; ivar < nVar; ivar++)
    {
        if(var0[ivar] != var0[ivar] || lmd0 != lmd0)
        {
            cout << "\n !!! Error: SU2 has diverged. Now exiting... !!! \n" << endl;
            cout << "------------------------CGKSSolver::calcFlux2D---------------------------\n" << endl;
            cout << "var0 has nan value! \n" << endl;

            cout << "CG:"<< endl;
            for(int ii = 0; ii < 2; ii++)
            {
                cout << aCG[ii] << "  " << sfCG[ii] << "   "<< bCG[ii] << endl;
            }

            cout << "tCG:"<< endl;
            for(int ii = 0; ii < 2; ii++)
            {
                cout << nodeCG[0][ii] << "  " << surfaceCG[ii] << "  "<< nodeCG[1][ii] << endl;
            }

            cout << "normal" << endl;
            cout << aNormal[0] << "  " << aNormal[1] << endl;

            cout << "aVar, bVar:"<< endl;
            for(int iivar = 0; iivar < nVar; iivar++)
            {
                cout << aVar[iivar] << "  "<< bVar[iivar] << endl;
            }

            cout << "nodeVar0:"<< endl;
            for(int iivar = 0; iivar < nVar; iivar++)
            {
                cout << nodeVar[0][iivar] << "  "<< nodeVar[1][iivar] << endl;
            }
            cout << " Gradient" << endl;
            for(int iivar = 0; iivar < nVar; iivar++)
            {
                cout << GradientA[iivar][0] << "  " << GradientA[iivar][1] << "  "
                    << GradientB[iivar][0] << "  " << GradientB[iivar][1] << endl;
            }

            cout << "var0:"<< endl;
            for(int iivar = 0; iivar < nVar; iivar++)
            {
                cout << surfaceVar[0][iivar] << "  "<< var0[iivar] << "  " << surfaceVar[1][iivar] << endl;
            }

            cout << "tau = " << tau  << "  deltaT = " << deltaT << endl;
            cout << "mu0 = " << mu0<< " eddyVis = " << eddyVis<< endl;
            exit(EXIT_FAILURE);

        }
    }


    //***********************************************************************************************
    // calculating the aBar[iParent]
    // aBar[0]: aBar l x
    // aBar[1]: aBar r x
    // aBar[2]: aBar y

    for(int iParent = 0; iParent < nParent; iParent++)
    {
        // calculating the aBar[indexChild][ivar]
        double tempdx = surfaceCG[0] - nodeCG[iParent][0];
        for(int ivar = 0; ivar < nVar; ivar++) B[ivar] = (var0[ivar] - nodeVar[iParent][ivar]) * rRho0 / tempdx;
        solveEquation(aPrandtl_Lam, var0[0], var0[1] * rRho0, var0[2] * rRho0, lmd0, B, aK, aBar[iParent]);
    }

    calcUmVnXPsi(UPsi[0], 0,  VPsi[0], 0,  a[0][1], Ksi2[0], Ksi4[0], B);
    calcUmVnXPsi(UPsi[1], 0,  VPsi[1], 0,  a[1][1], Ksi2[1], Ksi4[1], tempB);
    for(size_t ivar = 0; ivar < 4; ivar++) B[ivar] = rRho0 * (B[ivar] * surfaceVar[0][0] + tempB[ivar] * surfaceVar[1][0]);
    solveEquation(aPrandtl_Lam, var0[0], var0[1] * rRho0, var0[2] * rRho0, lmd0, B, aK, aBar[2]);
    //for(size_t ivar = 0; ivar < 4; ivar++) aBar[2][ivar] = 0.0;
    //***********************************************************************************************
    // calculating the ABar
    // the coefficient
    double gamma0 = deltaT - tau * (1.0 - exp(-deltaT / tau));
    double gamma2 = (-deltaT + 2.0 * tau * (1.0 - exp(-deltaT / tau)) - deltaT * exp(-deltaT / tau)) / gamma0;
    double gamma4 = (deltaT * exp(-deltaT / tau) - tau * (1.0 - exp(-deltaT / tau))) / gamma0;
    double gamma5 = - tau * (1.0 - exp(-deltaT / tau)) / gamma0;

    //for(int ivar = 0; ivar < 4; ivar++) B[ivar] = 0.0;
    // the macro-variables of interface
    //
    // gamma2 * u * aBar
    calcUmVnXPsi(U0Psi[1], 1,  V0Psi, 0,  aBar[0], Ksi02, Ksi04, B);
    calcUmVnXPsi(U0Psi[2], 1,  V0Psi, 0,  aBar[1], Ksi02, Ksi04, tempB);
    for(int ivar = 0; ivar < 4; ivar++) B[ivar] = (B[ivar] + tempB[ivar]);
    calcUmVnXPsi(U0Psi[0], 0,  V0Psi, 1,  aBar[2], Ksi02, Ksi04, tempB);

    for(int ivar = 0; ivar < 4; ivar++) B[ivar] = (B[ivar] + tempB[ivar]) * gamma2 * var0[0];

    for(int iParent = 0; iParent < 2; iParent++)
    {
        double aRho = surfaceVar[iParent][0];

        // (gamma4 + gamma5) * u * a
        calcUmVnXPsi(UPsi[iParent], 1,  VPsi[iParent], 0,  a[iParent][0], Ksi2[iParent], Ksi4[iParent], tempB);
        double tempB1[4];
        calcUmVnXPsi(UPsi[iParent], 0,  VPsi[iParent], 1,  a[iParent][1], Ksi2[iParent], Ksi4[iParent], tempB1);
        for(int ivar = 0; ivar < 4; ivar++) B[ivar] += (tempB[ivar] + tempB1[ivar]) * (gamma4 + gamma5) * aRho;

        // gamma5 * A
        calcUmVnXPsi(UPsi[iParent], 0,  VPsi[iParent], 0,  A[iParent], Ksi2[iParent], Ksi4[iParent], tempB);
        for(int ivar = 0; ivar < 4; ivar++) B[ivar] += tempB[ivar] * gamma5 * aRho;
    }


    for(int ivar = 0; ivar < 4; ivar++) B[ivar] *= rRho0;
    solveEquation(aPrandtl_Lam, var0[0], var0[1] * rRho0, var0[2] * rRho0, lmd0, B, aK, ABar);
    /*

       cout << "aBar[0]" << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[0][0] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[0][1] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[0][2] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[0][3] << endl;

       cout << "aBar[1]" << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[1][0] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[1][1] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[1][2] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[1][3] << endl;

       cout << "aBar[2]" << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[2][0] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[2][1] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[2][2] << endl;
       cout << setiosflags(ios::scientific) << setprecision(20)<< aBar[2][3] << endl;

       cout << "ABar" << endl;
       cout << setiosflags(ios::scientific) << ABar[0] << endl;
       cout << setiosflags(ios::scientific) << ABar[1] << endl;
       cout << setiosflags(ios::scientific) << ABar[2] << endl;
       cout << setiosflags(ios::scientific) << ABar[3] << endl;
       cout << endl;
       */

    //*********************************************************************************************************
    // flux
    // the macro-variables of interface

    // initializing the flux
    double* flux = new double[4];
    double* tempf = new double[4];

    //// g0
    tempf[0] = U0Psi[0][1];
    tempf[1] = U0Psi[0][2];
    tempf[2] = U0Psi[0][1] * V0Psi[1];
    tempf[3] = 0.5 * (U0Psi[0][3] + U0Psi[0][1] * V0Psi[2] + U0Psi[0][1] * Ksi02);

    double t00 = -tau * (exp(-deltaT / tau) - 1.0);
    double t01 = -tau * (deltaT * exp(-deltaT / tau) - t00);
    double t1 = deltaT - t00;
    for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] = tempf[ivar] * t1 * var0[0];

    double U0 = var0[1] * rRho0;
    double V0 = var0[2] * rRho0;
    // aBar * g0
    double t2 = -tau * t1 + t01;
    double tempf1[4], tempf2[4];
    calcUmVnXPsi(U0Psi[1], 2,  V0Psi, 0,  aBar[0], Ksi02, Ksi04, tempf);
    qs += tempf[nVar-1];
    qs -= U0*tempf[1];
    qs -= V0*tempf[2];

    calcUmVnXPsi(U0Psi[2], 2,  V0Psi, 0,  aBar[1], Ksi02, Ksi04, tempf1);
    qs += tempf1[nVar-1];
    qs -= U0*tempf1[1];
    qs -= V0*tempf1[2];

    calcUmVnXPsi(U0Psi[0], 1,  V0Psi, 1,  aBar[2], Ksi02, Ksi04, tempf2);
    qs += tempf2[nVar-1];
    qs -= U0*tempf2[1];
    qs -= V0*tempf2[2];

    for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] += (tempf[ivar] + tempf1[ivar] + tempf2[ivar]) * t2 * var0[0];

    // ABar * g0
    calcUmVnXPsi(U0Psi[0], 1,  V0Psi, 0,  ABar, Ksi02, Ksi04, tempf);
    qs += tempf[nVar-1];
    qs -= U0*tempf[1];
    qs -= V0*tempf[2];

    double t3 = 0.5 * deltaT * deltaT - tau * t1;
    for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] += tempf[ivar] * t3 * var0[0];

    // g / a
    for(int iParent = 0; iParent < 2; iParent++)
    {
        double aRho = surfaceVar[iParent][0];

        // g
        //
        tempf[0] = UPsi[iParent][1];
        tempf[1] = UPsi[iParent][2];
        tempf[2] = UPsi[iParent][1] * VPsi[iParent][1];
        tempf[3] = 0.5 * (UPsi[iParent][3] + UPsi[iParent][1] * VPsi[iParent][2] + UPsi[iParent][1] * Ksi2[iParent]);

        for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] += tempf[ivar] * t00 * aRho;

        // u * a * g 
        calcUmVnXPsi(UPsi[iParent], 2,  VPsi[iParent], 0,  a[iParent][0], Ksi2[iParent], Ksi4[iParent], tempf);
        calcUmVnXPsi(UPsi[iParent], 1,  VPsi[iParent], 1,  a[iParent][1], Ksi2[iParent], Ksi4[iParent], tempf1);
        for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] += -(tempf[ivar] + tempf1[ivar]) * (t01 + t00 * tau) * aRho;
    }

    // A * g * t00 
    for(int iParent = 0; iParent < 2; iParent++)
    {
        double aRho = surfaceVar[iParent][0];

        // A * g * t00
        calcUmVnXPsi(UPsi[iParent], 1,  VPsi[iParent], 0,  A[iParent], Ksi2[iParent], Ksi4[iParent], tempf);

        for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] += -tempf[ivar] * tau * t00 * aRho;
    }

    // Prandtl number fix
    qs *= -tau * deltaT * var0[0];

    //flux[nVar-1] += (1.0 / aPrandtl_Lam - 1.0) * qs;

    for(int ivar = 0; ivar < nVar; ivar++) flux[ivar] *= surfaceVolume;

    // transform the coordinate system
    //
    // flux
    temp1 = flux[1]; temp2 = flux[2];
    flux[1] = temp1 * aNormal[0] - temp2 * aNormal[1];
    flux[2] = temp1 * aNormal[1] + temp2 * aNormal[0];

    delete[] tempf;

    for(int iParent = 0; iParent < nParent; iParent++)
    {
        for(int ivar = 0; ivar < nVar; ivar++)
        {
            delete[] GradientArr[iParent][ivar];

        }

        delete[] GradientArr[iParent];
    }

    delete[] GradientArr;

    return flux;
}

void solveEquation(double aPrandtl, double aRho, double aU, double aV, double aLambda, double* B, double aK, double* aX)
{
    double R4 = 2.0 * B[3] 
        - (aU*aU + aV*aV + 0.5*(aK+2.0)/aLambda) * B[0];
    double R3 = B[2] - aV * B[0];
    double R2 = B[1] - aU * B[0];

    aX[3] = 4*aLambda*aLambda * (R4 - 2 * aU * R2 - 2*aV * R3);
    aX[3] /= (aK+2.0);
    aX[3] /= aPrandtl;
    aX[2] = 2 * aLambda * R3 - aV * aX[3];
    aX[1] = 2 * aLambda * R2 - aU * aX[3];
    aX[0] = B[0] - aU*aX[1] - aV*aX[2] 
        - 0.5*aX[3]*(aU*aU + aV*aV + 0.5*(aK+2.0)/aLambda);
} //: solveEquation


void transformCoordSystem(double* aCG, double* aVar, double** aGradient, 
        double* rCG, double* rVar, double** rGradient, double* normal, double* tangent)
{
    transformCoord(aCG, rCG, normal, tangent);
    transformConsVar(aVar, rVar, normal, tangent);
    transformConsGrad(aGradient, rGradient, normal, tangent);
}

void transformConsVar(double* aVar, double* rVar, double* normal, double* tangent)
{
    rVar[0] = aVar[0];
    rVar[3] = aVar[3];
    rVar[1] = aVar[1] * normal[0] + aVar[2] * normal[1];
    rVar[2] = aVar[1] * tangent[0] + aVar[2] * tangent[1];
}

void transformConsGrad(double** aGradient, double** rGradient, double* normal, double* tangent)
{
    double temp1 = aGradient[0][0];
    double temp2 = aGradient[0][1];
    rGradient[0][0] = temp1 * normal[0] + temp2 * normal[1];
    rGradient[0][1] = temp1 * tangent[0] + temp2 * tangent[1];

    temp1 = aGradient[3][0];
    temp2 = aGradient[3][1];
    rGradient[3][0] = temp1 * normal[0] + temp2 * normal[1];
    rGradient[3][1] = temp1 * tangent[0] + temp2 * tangent[1];

    temp1 = aGradient[1][0];
    temp2 = aGradient[1][1];
    double temp3 = aGradient[2][0];
    double temp4 = aGradient[2][1];

    rGradient[1][0] = normal[0] * normal[0] * temp1 
        + normal[0] * normal[1] * temp2
        + normal[0] * normal[1] * temp3
        + normal[1] * normal[1] * temp4;

    rGradient[2][0] = normal[0] * tangent[0] * temp1 
        + normal[1] * tangent[0] * temp2
        + normal[0] * tangent[1] * temp3
        + normal[1] * tangent[1] * temp4;

    rGradient[1][1] = tangent[0] * normal[0] * temp1 
        + tangent[1] * normal[0] * temp2
        + tangent[0] * normal[1] * temp3
        + tangent[1] * normal[1] * temp4;

    rGradient[2][1] = tangent[0] * tangent[0] * temp1 
        + tangent[0] * tangent[1] * temp2
        + tangent[0] * tangent[1] * temp3
        + tangent[1] * tangent[1] * temp4;

}

void transformCoord(double* aCoord, double* rCoord, double* normal, double* tangent)
{
    rCoord[0] = aCoord[0] * normal[0] + aCoord[1] * normal[1];
    rCoord[1] = aCoord[0] * tangent[0] + aCoord[1] * tangent[1];
}



